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We study a model, motivated by a bioremediation process, describing a cross-diffusion movement of a 
bacteria population b attracted by a chemoattractant signal c, in a nonhomogeneous stratified medium with 
n layers. We assume that this reaction-diffusion process is characterized by a low rate of degradation and 
a low diffusion coefficient of the chemoattractant, expressed in the model by a small parameter s. The 
model consists of n systems of nonlinear parabolic equations with transmission conditions between layers. 
We prove a global-in-time solution for the asymptotic model setup with respect to the small parameter of 
the problem, for arbitrarily large initial data. Next, we deal with the control problem focusing mainly on 
the reduction of the chemoattractant concentration, by acting upon the initial distribution of the bacteria 
population bo- To this end, we prove the existence of a solution to the control problem and determine the 
optimality conditions. 
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1. Introduction 

Chemotaxis is the biological process of directed movement of cells in response to a chemical signal 
emitted by a substance or by another population, called chemoattractant, and play an important role 
in the interaction of cells with their environment. Chemotaxis is a complex process which involves 
many aspects such as various species of cell populations b can interact; the movement can be either 
towards the higher concentration of the signal (positive chemotaxis) or away from it (negative 
chemotaxis); the chemical signal can be secreted by a population b itself, and not necessarily by 
an external source; the chemotactic process can lead to an aggregation of the attracted individuals 
and to a production of the chemoattractant, or contrary, it can determine the degradation of the 
chemoattractant. Consequently, the theoretical understanding of these processes determined a 
growing interest in their mathematical modelling, which at its turn raised challenging problems. 
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The origin of the fundamental model is given in the work of Patlak [24]. Later, Keller and Segel 
[ 1 8-20] introduced a similar model based on another assumptions. Since then, a rich mathematical 
literature on various versions of the model has been emerged, mainly focusing on the well- 
posedness of it and we refer the reader to a very comprehensive survey in [15]. 

In this paper, we shall denote the density of the cell population by b and the density of the 
population spreading the signal by c, both of them being functions of time t , and space variable 
f. 

A chemotactic system with only one population of cells consists of two equations for b and c 
with initial and boundary conditions: 



- V • (D(b,c)Vb) + V • (K(b,c)bVc) 


= g(b,c)-h(b,c), inQ, 


(1) 


b(0,^) 


= b 0 (%), in £2, 


(2) 








^-V-(5(Z>,c)Vc) 
at 


= <p\(b,c) - (p 2 {b,c), inQ, 


(3) 


c(0,f) 


— cq, in £2. 


(4) 



Here, the time runs in (0, T), T is finite, f e Q which is an open bounded subset of R d and 
Q — (0,T) x Z2. In these equations, D(b,c) and S(b,c) represent the diffusion coefficients of 
the attracted population b and chemoattractant c, respectively; gib, c) and h(b, c) are functions 
describing the rates of growth and death of b; and q>\ (b, c) and (p2(b, c) stand for the production 
and degradation of the chemoattractant. Further, we denote the right-hand sides by f(b,c) — 
g(b, c) — h(b, c) and <p(b, c) — q)\ (b, c) — (p 2 (b, c). They may depend on the space variable, too. 
The function K characterizes the chemotactic sensitivity and the cross-diffusion term in the first 
equation is indeed able to enforce the spontaneous emergence of structures provided that the 
process of chemotactic migration is accompanied by a production of the signal substance by the 
cells themselves [14]. Thus, the cross-diffusion term and the kinetic term <p(b,c) in the second 
equation can lead to the blow-up of the solution, even if the growth-death rate f(b, c) is zero and 
D, S, K are constant. In the literature, the chemotactic system has been approached in simplified 
versions, able to avoid blow-up and to allow global solutions. We cite some more recent results. 

In the l-D case, it has been shown that blow-up does not occur [16] for D — 1, K constant 
and S — s small. When the space dimension d is greater or equal to two, the solutions generally 
exhibit blow-up, this being influenced by the model parameters and the characteristics of the initial 
data [13,14,25]. For example, in [17], a chemotaxis motion with constant diffusion coefficients 
is studied by using a nonlocal gradient sensing term to model the effective sampling radius of 
the species. In [9], Dyson et al. use a nonlocal term to model the species-induced production 
of the chemoattractant, (p 2 (b,c,^), in order to prevent blow-up in the d-D space, considering 
that the diffusion coefficients are constant. They prove the existence of solutions, which exist 
globally, and are L°°-bounded on finite time intervals. In [28], the system (l)-(4) is considered with 
D(b,c) — S(b,c) — K(b,c) — l,f(b,c) — 0and^(Z?,c) = —foe, and it is shown that for arbitrarily 
large initial data, this problem admits at least one global weak solution for which there exists 
T > 0, such that (b, c) is bounded and smooth in (T, oo) x £1. The paper [8] deals with two types of 
reaction-diffusion systems, one arising in chemotaxis and the other in angiogenesis. The first refers 
to Equations (l)-(4) with D(b, c), S(b, c) and K(b, c) constant,/(fo, c) = 0, cp(b, c) — b — ac and 
the equation for c stationary, that is, — Ac — b — ac. This equation is obtained as a case limit, for 
e — > 0, of Equation (3) in which cp\ (b, c) — (p 2 {b, c) — ( l/e)(b — c), and the diffusion coefficient 
S is of the order of 1/e, with £ small. For this system (as well as for that of angiogenesis type), it 
is shown in [8] that when the L d ^ 2 norm of initial data is small enough (for d > 2), then there is a 
global (in time) weak solution that stays in all LP spaces with max{l;cf/2 — 1} < p < oo. In [4], 
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the same system but with a — 0 is studied in R 2 and a detailed proof of the existence of weak 
solutions below the critical mass, above which any solution blows up in finite time in the whole 
Euclidean space, is given [26]. The stability of the stationary solutions to a chemotaxis system 
was proved in [1 1] for D — 1, S(b, c) — 0 and a general function cp(b, c). In [23], existence and 
uniqueness are studied for slow and singular fast diffusion of the cells in the case with a stationary 
equation for the chemoattractant. 

Besides applications in biology, a chemotaxis model can apply in environment bioremediation 
in the cases when specific bacteria b are injected into a polluted medium (soil or water) with the 
purpose of cleaning it from an inside spread pollutant c [5,10,27]. 

Our study is motivated by an application to environment bioremediation and focuses on the 
case in which the kinetic term and the diffusion coefficient of the chemoattractant (pollutant) c 
have a weak influence on the flow, meaning that the rate of degradation of the chemoattractant is 
slow and it diffuses very little (or not at all, as in the case of oil polluting an environment). 

Roughly speaking we shall start from a model reading as 

d ± _ V • (D&Vb) + V • (K(b,c,i;)bVc) =f(b,c,$), in Q, 
at 

b(0,!;) = b 0 ($), in £3, 

^-sV -(8(i;)Vc)=ecp(b,c,i;), mQ, 
at 

c(0,§) =c 0 (§), inC2, 



where e is a small parameter in front of the diffusive and the kinetic terms for the chemoattractant. 
Such a model is obtained by making dimensionless Equations (l)-(4). Moreover, we assume that 
at the initial time t — 0 the chemoattractant concentration c 0 is constant and that there exists a 
certain nonlinear dynamics growth-death of the cell population, indicated by/. 

The model will be implemented into a nonhomogeneous medium, fact evidenced later by 
a particular space dependence of the coefficients. The aim is to investigate the way in which 
the solution depends on the small parameter e and not to study the limit model when s — > 0. 
Accordingly, we shall not pass to the limit, but use a perturbation technique [7], by which the 
solution is expanded in series with respect to the powers of the small parameter, and retain the 
systems of s m -order of approximation, obtained by equating the coefficients of e m in Equations (5). 

The main goal of the paper is to study the possibility of controlling the environment cleaning 
by acting upon the initial distribution of the bacteria. More exactly, the intention is to design 
the initial distribution of the bacteria such that the chemoattractant mean concentration should 
decrease under a certain critical value c crt within a time period T. 

A secondary task of the control problem is motivated by the fact that bioremediation can be 
problematic because sometimes it might become difficult to remove the microorganisms released 
into the environment in order to let them at acceptable levels in the soil [21]. Consequently, 
for enhancing a complete and efficient decontamination process and for avoiding higher costs 
of a further cleaning of the environment from the remained bacteria, we investigate by the same 
control problem whether a moderate bacteria proliferation, whose growth is limited by a prescribed 
density value b c , would be sufficient for achieving the main objective of reducing the pollutant 
concentration up to a value under the dangerous critical threshold c crt . To this end, the state system 
is studied in the framework of an appropriate functional setting and a global-in-time solution for 
the asymptotic system derived from Equations (5) is obtained. Then, the existence of at least a 
solution to the control problem is proved and the optimality conditions proving information about 
the initial distribution of the bacteria bo are computed. 
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2. Statement of the model and the perturbation technique 

We consider a nonhomogeneous 3D right cylinder domain 

£2 = = (x,y,z) el 3 ;i€ (x 0 ,x n ), f = (y,z) e 

where the base £2' is an open bounded subset of R 2 , with a boundary of class C 2 . The medium 
nonhomogeneity is modelled by a stratification of £2 in n parallel layers along the Ox axis, the 
separation of the layers being determined by the different values that movement parameters may 
have in each layer. More exactly, a stratification can be put into evidence when certain parameters 
do not depend on x (here x is the stratification variable) or are constant in each layer (x,_ i , x,) but 
have different values from a layer to the other. 

Therefore, the domain £2 consists of n subdomains Qi, having the boundaries 3 £2, = r,_i U 
r ; U r] at , i — 1, . . . , n, where r! at are the lateral boundary of £2,- and T ; = {x — x ; }, i = 0, . . . , n. 
The surfaces To and V n are the external horizontal boundaries, while T, with i — 1, . . . , n — 1 are 
the boundaries between layers. We denote 

Qi := (0, T) x S2,-, Ej := (0, T) x r { , £, lat := (0, T) x r, lat , i = 1, . . . , n. 

In each layer the chemotaxis process is modelled by two equations, one for b t and the other 
for the chemoattractant c,. The interaction between the layers is established by transmission 
conditions for that is, the continuity of the solutions and fluxes. We assume that the system is 
closed for bi and c,, namely the fluxes across the exterior frontiers are zero. The chemoattractant 
c; may display jumps at the interface between layers, because the ^"-approximation systems will 
not be anymore of diffusion type, as we shall see. Therefore, there is no need to specify these 
conditions. 

The stratification is set by the values for £),-, <5 ; , c,- q assumed constant in each layer i, but different 
for two consecutive layers. Also, the expressions of the functions/;', Ki, (p i7 which are assumed 
not to depend explicitly on §, are different from one layer to another. 

With these considerations, the dimensionless mathematical model for a chemotaxis movement 
in a nonhomogeneous stratified medium has n equations for the unknowns bi and a (the density of 
the cell population and the chemoattractant concentration, respectively, in each layer i — 1, . . . , n) 
and reads [1] 

^ - DDjAbi + KV ■ [bjKiibi, cOVc,] =ffi(b h c,-), (f, I) e Q h (6) 
at 

dci 

— - SSiAQ = q> Vi {bi, a), (t, H) e Qi, (7) 
of 

c i (0,|) = c M) , feO,-, (8) 
i i (0,|) = ft ii0 (|), fe£2i, (9) 

for all (' = 1, . . . , m, where c,,o (assumed constant) and fe,-,o(t) are initial conditions for c, and fo,. 
At the interface between two layers, we have the conditions 

obi - 3c,- - db i+l - dci + i 

-DDi-^+KbiKiib^a)-^ = -DDi +x ^ + Kbi +l K i+x {b i+u Ci +l )^L on Ej , (10 ) 
dx 9x dx dx 

fo,=^+i onS ; , (11) 
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db, 
DDi — - 

dx 


+ Kb 1 K i (b u c l )— 1 - 
dx 


= 0 


on So, 






db„ 
DD n — - 

dx 


dc„ 

+ Kb n K n (b n ,c n )-^ 
ox 


= 0 


on E„, 








Vbi ■ v 


= 0 


on E, lat , 


i=l,.. 


• ,n, 




dc\ 
dx 


= 0 


on So, 








dc n 
dx 


= 0 


on E„, 








Vet ■ v 


= 0 


on £, lat , 


i=l,.. 


. ,n. 



for i — 1, . . . , n — 1, and on the exterior horizontal and lateral boundaries, we set 

(12) 

(13) 
(14) 

(15) 

(16) 
(17) 

Here, v is the unit outer normal to r! at and d/dv is the normal derivative. 

Generally, c, and its normal derivative 8j(dci/dx) can have a jump when crossing the internal 
boundaries T,. As we said, do not specify here these conditions because they will not intervene 
in the asymptotic model we shall study. 

We stress that this model is presented in a dimensionless form and the constants D, K,f, 8, <p are 
dimensionless parameters. They appear after the dimensionless procedure by which all dimen- 
sional variables and function are scaled by certain characteristic values indicated by the subscript 
a (e.g. L a for length, T a for time, D a for diffusion coefficient, etc.). We do not present the dimen- 
sionless computation, but we assert that by this procedure the dimensionless parameters are given 
by 

r. T a - T a - c a T a - T a _ T a 

D = -^D a , S = — 2 b a , K = —^Ka, f = -f a , <p = —<p a . 

"a c a 

If by such a procedure the computed values 8 and <p follow to have the same order of magnitude, 
much smaller than that of the other parameters, we set in Equations (6)— (17) 

(p = e,8 = s (18) 

and consider the other dimensional parameters of order 0(1) as e — > 0. We agree that this model 
characterizes the chemotactic process with a chemoattractant (e.g. oil in a bioremediation process) 
that diffuses very slow, and it is degraded with difficulty by the bacteria. This process happens in 
a different way in each layer, due to the nonhomogeneity modelled by the various functions/;-, 
Kj, (pi, the structure of the initial data and different diffusion coefficients. 



2.1. Hypotheses 

To approach system (6)— (17) we assume the following hypotheses for all i — 1, . . . , n : 

(ii) Ci t o is constant, c, o > 0, and there exist an i such that c, o > 0; 
(h) bifi > 0, and there exist an i such that bi$ > 0; 

(1 3 ) A > D 0 > 0, D 0 = min, =I ^ A; 

(14) (^1,^2) — > Ki(r\,r2) are of class C 1 with respect to r\, \Ki(r\,r2)\, |(3^,/3ri)(ri, r{)\ are 
bounded, Vri, r-i e R; 

(is) (r\,r-i) — > (fi(r\,r2) are of class C 2 with respect to n and \<Pi(r\,r2)\, \(d<Pi/dr{)(ri, r?)\, 
\(d 2 (pi/drl)(ri,r2)\ are bounded. 

We remark that equations with nonlinear terms do not generally admit global solutions in 
time [12], but by the perturbation technique we shall deduce a global solution for this asymptotic 
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model, even under the assumption of a polynomial form for ft. Thus, for all i — 1, . . 
assume 



/i g C'(i x 1) Vri,r 2 eR, 
fi(rur 2 )n <0 Vn.fi eR, 
./■(0,r 2 )=0 Vr 2 eM, 

3/- 



9r 2 
8ri 



(n,r 2 ) 



<Cj(l + |rin Vn.neR, 



with CJ, k' c positive for i = !,...,«, where 



0</?<2 if*/ = 3, 
0</?<oo if = 1,2, 



,«, we 

(19) 
(20) 
(21) 

(22) 
(23) 



and d is the space dimension. 



2.2. e a -order and s 1 -order approximations 

We write the series expansions of all functions, with respect to the small parameter e = cp — S. 
Denoting generically </>,(f, §) and $;(&,■, c,), we have 



<Pi(t,i;) = <P?(t,I;) + e<f>l(t,$) 



*«(*i,c,-) = *,-(*?,c?) + e(*i)*»(*?,c?)* I ? + c?)c? + ■ ■ ■ 



(24) 



where 0™ corresponds to the s m -approximation of 0, (m = 0, 1, . . .), and (4> ( )fe,, (3>s)c ( represent 
the derivatives of 3>, with respect to bi and c, . 

We replace these series in the system (6)— (17) and by equating the coefficients of the mm 
powers of s, we deduce the systems corresponding to e m -order approximation. 

The e°-order approximation for a reads 



9c° 

c?(0,f) = c w inn,-, 
for each i = implying that c° is a positive constant in each layer, 



For we get 



8b[ 
dt 



q(t, |) = c,-,o = constant V(f, § ) G g ; , *'=!,...,«. 



-DDiAb* =ffi(b-,Ci,o) in ft, 

&?(0,£)=&y>(£) in ^> i=l,...,n, 



(25) 
(26) 

(27) 



(28) 
(29) 
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DD,-^(a) = Dfl i+l ^(a) ; (»,f)eSi, i=l,...,B-l, (30) 
dx dx 

b$(.t,ft = b? +l (f,S), (M)eS,-, i = l,...,n-l, (31) 

-r^(M) = 0, (M)eS 0 . (32) 

-^(f,t) = 0, (a)eE„, (33) 
dx 

Vft°-v = 0, (a)eE, lat , i=l,...,n. (34) 

Further, identifying the coefficients of e 1 and using Equation (27), we obtain the system for the 
e 1 -order approximation, 



<Pi(b°, c ( -,o) in Q h /=!,..., n, (35) 



3c? 
3t 

c?(0,£)=0 in £2,-, /=l,...,n, (36) 

3ft 1 

-f- - DDjAb] -f(fi) bi ,(bl Cija)b] = F f (f, $) in g„ (37) 
or 

ft 1 (0,|) =0 in £2,-, i= l,...,n, (38) 



/ . 9ft, 1 - 3ftL,\ 

= Gi(t,Xi,?) on £,-, 



(39) 



ft 1 = fti +1 on E„ ; = 1, . . . , n - 1, (40) 

—DD\^p- = Go(t,Xo,%') onEo, (41) 
3x 

3ft 1 

DDn— 1 — G„(t,x n ,%') onE„, (42) 
3x 

Vft 1 ■ v = 0 on £ lat , i = 1, . . . , n, (43) 
where we have denoted £' := (y,z) e £2', 

= (ftfdctib^c^cl - KV ■ [ft°(f,?)^(^,Q,o)Vc 1 ])(f,?), (44) 

for ;' = 1, . . . ,«, 

Gi(t,Xi,%') = K I -^.(ftO )C . 0 )_L + ftO +i ^. +l( fcO +i)C . +lo) _L±! j (45) 
for z = 1 , . . . , n — 1 , 

G 0 (f,x 0 ,r) = * (ft°, Cl ,o)^ (*,*>,§'). (46) 

G n {t,x n ,^') = K (blK n {bl,c n $) d -^j (47) 

After solving the system for the e°-order approximation, c] follows by Equations (35) and (36) 
and so the functions F,(f, f ), Gi(t,Xi, £' ) are known. 
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The next approximations (for m > 2) lead to systems having similar forms as that for the e - 
order approximation, and so they pose the same mathematical problems. That is why we do no 
longer write them. 

Definition 2. 1 We call an asymptotic solution to Equations (6)-(17), up to the order of 
approximation s 2 , a pair of functions (bj, c ; ), i — 1, . . . , n, 

b i {t,H)^b () i (t^) + £b)(t^), 

(48) 

~C i {t,t;) = c a i (t,$;) + Sc](t,H), 

where c° and c] are the solutions to the differential equations (25) and (26) and (35) and (36), 
respectively, and b° t and b) are the solutions in the sense of distributions to Equations (28)-(34) 
and (37)-(43). 

We shall rigorously explain the definition of the solutions and b) in the next sections, in 
which the well-posedness of these systems will be studied. 

2.3. The control problem 

As we have already explained, the aim of the control problem is to provide information about 
the necessary initial density of bacteria, bo(%) = (&j,o(t))s=i,...,n> such that two objectives would 
be achieved. The main one is related to the environment cleaning that is, to force the decrease 
of the chemoattractant mean concentration under a critical threshold c crt . This means to mini- 
mize in the cost functional the mean positive part of the difference between the concentration 
c(?,|) — c°(f, £) + scj (t, f ) given by Equation (48) and c crt . The second objective aims at realiz- 
ing the necessary pollutant concentration decrease by a process limiting a too large proliferation 
of the bacteria. Consequently, in the cost functional we add a term expressing the restriction of 
the bacteria growth, by minimizing the mean positive part of the difference between b®(t,%) and 
a constant prescribed value b c . We introduce this control problem to give a quick response (rather 
than an extremely accurate one) and to indicate a first decision for a process evolving in real time. 
That is why we accept that the consideration of the s°-order approximation b ( - only, which is the 
dominant term in the asymptotic expansion of £>,, is motivated. A more accurate computation, if 
necessary, may take into account the further approximations. 
Mathematically, we have to minimize the cost functional 

i " r 

J(b lfi , b hn ) =-J2 ((C°(f, % ) + EC\ (t, f ) - C crt )+) 2 d£ df 

z ,=i J Q< 

n . 

+ ^J2 a^(t^)-b c ) + ) 2 ^dt, 

where the superscript '+' means the positive part (i.e. <p + (i; ) = max{<p(£), 0}). A certain choice 
of the constant a can induce a lower (or a higher) influence of the first term with respect to the 
second one in the cost functional. In particular, if the limiting condition of the proliferation of the 
bacteria is disregarded, the constant a can be set zero. 
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Therefore, the control problem is 

Min J{b h0 ,...M,n) CP) 

{b w eU,} i=1 „ 

subject to Equations (28)-(34) and (35)-(36), where 

U, = {b U o e L~(£2,-); hm < bud!) < b m }. 

We recall that in our model, the initial distribution of the density of the bacteria may depend on the 
space variable § (i.e. &o(£) = (fr;,o(£));=i,...,n) an d omv the chemoattractant initial concentration 
is considered constant (cq = (c, o)i=i,. ..,«)• A natural restriction for the initial densities of the 
bacteria b ;> o (i — 1, . . . , «), representing the controllers in the control problem, is to impose their 
boundedness. The minimum and maximum bounds, Z?,; m and Z?, iM , are fixed numbers, 0 < b^ m < 
bt y M- However, as we see later, this restriction does not necessarily imply that the minimum of / 
must be reached at fc, m . 

Finally, to be consistent with the requirement expressed by the first term in J, the maximum of 
the initial datum is chosen less than or equal to b s , the appropriate upper threshold envisaged for 
the solution b ( -. In the case with a — 0, the choice of Z?, iM is free. 



3. Well-posedness of the state system 

In order to approach the control problem we need some information about the e° and e Estate 
systems. We begin by a few preliminaries. 



3.1. Preliminaries and functional setting 

Since by our assumptions, c, o remains constant in each layer, one notes that in Equation (28) 
depends only on b® and so one can denote 

Hi(T) := -Mr, c u o) VreK, i = 1, . . . , n. (49) 

The assumptions (19)— (23) imply that the functions /Lt, e C^M), 

0 < fXi(r)r Vr e K, (50) 
/*,-(0)=0, (51) 
\Hi(r)\<fC[(l + \rn VreM, (52) 

with p G [0, 2) if d = 3 and p > 0 if d = 1 or d = 2. 

We rewrite the system (28)-(34) for the e°-order approximation without indicating the 
superscript '0' 

^i-DD i Ab i + p i (b i )=0 in ft, (53) 
at 

i i (0,t) = i ii0 «) in^„ i=l,...,n, (54) 

dbi dbi+i 
A-^=A+i^ ±i onE„ (55) 
ox dx 
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b t — b i+ \ on E„ i = 1, . . . , n - 1, 

3x 
3^ 
3x 



= 0 on E 0 , 



0 on E„, 



V*, . v = 0 onEf, i = !,...,«. 



(56) 
(57) 

(58) 
(59) 



We introduce now the global functions denned in the following way (indicated by a generic 
notation <p,<po)'- 



0„(f,£), X € (x„_l,X n ), 
01,o(£), * G (jfO.JCl) 

0o(x) = : 

_0n,o(£)> * G (x n _uXn) 

(for a function which is constant in each layer, as cq or D), 



(60) 



3>(*,c,x) = 



4>i(fei,ci), x e (x 0 ,xi), 



$«(*«, c„), x e (x„_i,x„) 



(61) 



(for a function depending on b and c), where we recall that § = (x,y, z). We mention that we 
include the constant D in the definition of D(x). 

We note that assumptions (50)-(52) and (ii)-(is) imply similar properties for the functions 
defined before. Namely, we have 



0 < fi(r,x)r Vr e R, a.e. x e (xo,x„), 
\li' r (r,x)\ < Ci(l + \rf) VreM, 
/x(0,x) = 0 



(62) 
(63) 
(64) 



uniformly with respect to x, where C\ =/max (= i > and /x' r is the derivative with respect 
to r, 



co>0, 

hiH) > o, 

D(jt) > D 0 > 0. 



(65) 
(66) 
(67) 



We recall that we consider p e [0, 2) if J = 3 and /? > 0 if d — 1 or <i = 2. 

System (53)-(59) will be treated in the functional framework of the Sobolev space V — H l (£2) 
endowed with the standard norm H^lly = (||Vi 2 + ||Vi/i 2 ) 1/2 , and its dual V, with the pivot 
H = L 2 (f2), such that V C H c V . The value of g e V at f e V is 



git) = (g,t)v,v, 
where (•, -)y,v represents the duality between V and V. 



(68) 
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We specify that for the writing simplicity we shall denote the scalar product and the norm in 
i 2 (£^) by (-,-) and || - |l- 

We define the operator Aq : V — > V by 



(A 0 b, $)y,y = V / [DDiVbi ■ Vf + m{bi)f] d£ 

i=l Jn > 

= f[D(x)Vb- Vf + i±(b,x)f]dl; ViA e V (69) 

and restrict it to L 2 (Q) by the operator A : D(A) C L 2 (Q.) ->• L 2 (£2), 

A6=A 0 £ Vfce£»(A) (70) 

with D(A) = {be V,Ab e L 2 (Q)}. 

Thus, the functional abstract setting of our problem is 

db 

— (t)+Ab(t)=0 a.e.fG(0J), (71) 
at 

b(0) = b 0 . (72) 
3.2. Existence for the -order approximation 

The well-posedness for Equations (71) and (72) is concentrated in the following. 

Theorem 3.1 Assume conditions (62)-(64), and let bo e D(A). Then, problem (71) and (72) has 
a unique solution 

b e W 1 ' oo ([0, 71; L 2 (£2)) n L°°(0, 7 1 ; V) D C([0, 71; L 2 (£2)). (73) 
77ie solution satisfies the estimate 

IIMOIIv < C v , M«»yfe[0J], (74) 
where the constant Cy depends on the problem data. Moreover, 

bi G L 2 (0, T;H 2 (Q.i)), i=\,...,n. (75) 
Finally, let bo > 0 a.e. m £2. 7%en, f&e solution b to Equations (71) awd (72) satisfies 

b(t)>0 a.e.inQ Vt e [0,T]. (76) 

Proof The proof is split into two steps and is given in [1]. Further, we indicate the arguments. 
We assert by the properties (62)-(64) and by handling appropriate Sobolev inequalities that the 
operator b — > fi(b,x) turns out to be locally Lipschitz from V to L 2 (Q), uniformly in x. In the 
first step, according to the method presented in [3], we reduce the problem to a case with a globally 



Journal of Biological Dynamics 



99 



Lipschitz operator, by approximating n{-,x) by 

ix(b,x), 



fi N (b,x) - 



ix 



(—A 

\\\b\\v ) 



ll*(0llv<#. 
\\b(t)\\ v >N, 



(77) 



for each TV natural. This truncated operator is Lipschitz from V to L 2 (£2) for each N fixed, and 
we consider the approximating problem 



df 



■(f) +A N b N (t) = 0 a.e. t G (0, T), 
MO) = 6o, 



(78) 
(79) 



where A w is the operator A (defined by Equation (70)) in which fi(b,x) is replaced by fi N (b,x) 
and D(A N ) — D(A). 

For each N, it can be proved that A N is quasi-m-accretive on L 2 (£2) and then, for bo e D(A N ) — 
D(A), problem (78) and (79) has a unique strong solution 

b N G ff'^dO, T]; L 2 (Q)) n L°°(0, T; D(A)) fl C([0, 1]; L 2 (fi)). 

Some further necessary estimates are deduced for this solution, that is, 

II MO II < 11*0 II Vfe[0,r], (80) 

IIVMOII 2 < C 0 ||My, for any fe [0,71, (81) 



implying Equation (74), with Co depending on fixed data of the problem (p, T, Dq,Doo), and 

2 



db N 



At 



LHQ) 



<C 0 \\bo\\ 2 v . 



(82) 



The proof of the m-accretiveness of A N , as well as all estimates written before involve some 
technical computations given in detail in [1]. 

The positiveness is proved on the basis of Stampacchia's lemma and Gronwall's lemma and 
show that the solution falls within the accepted physical domain of positive densities. 

The regularity H 2 (a.e. t) of the restriction of b to each layer b- t comports a fine technique 
because it must adapt the known results [2,6] for a domain with a regular boundary 9 £2 to this 
transmission problem with n — 1 interfaces. These arguments are given in [1], too. 

Finally, Equation (74) opens the ways to the second step by going back to Equations (71) 
and (72). In fact, if we choose/? = Cy and N large enough, N > R, we get that A N b N (t) —AbN(t). 
Therefore, b N (t) with N large enough is the solution to problem (71) and (72), and further, it will 
be denoted by b. 

For proving the uniqueness of the solution, we consider two solutions b and b corresponding 
to the same initial data bo. By the previous proof if 

W> sup ||ft(f)||v+ sup \\b(t)\\v, 

te[0,T] te[0,T] 

then we get b(t) — b N (t) and b(t) — b N (t), where b N it) is the solution to Equation (78) and (79). 



It is obvious that in each layer the functions b t (meaning in fact the ^-approximation b®) have 
the regularity induced by b restricted to Q it for i — l,...,n, that is, 

b° G W 1 ' oo ([0,r];L 2 (^,))nL oo (0,r;// 1 (fi / ))nC([0,r];L 2 (^,))nL 2 (0,r;// 2 (^ / )). (83) 
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3.3. Existence for the e 1 -order approximation 

The functions c\ are computed by Equations (35) and (36), 

c\{t,$) = f W (*?(T,£),Cy,)dT, i=l,...,n, (84) 
Jo 

and by the boundedness hypothesis (i 5 ) and Equations (83) and (75), we derive that 

c ] e ff'-^fa) n ff'^flo, r];ff'(i2i)) n c'([0,r];L 2 (a)) n ff'^o, r];ff 2 (^i)). (85) 

Existence for the solution to Equations (37)-(43) for the e border approximation b]{t,%) is 
studied by the Lions' theorem for the time-dependent case [22] and the conclusions are [1] the 
following. 

Theorem 3.2 System (37)-(43) has a unique solution 

b 1 e W U2 ([0,T];V') n L 2 (0,T;V) n C([0,T];L 2 (tt)). (86) 



Corollary 3.3 Problem (6)— (17) admits a unique asymptotic solution up to the order of 
approximation s 2 , 

b e C([0, T];L 2 (£1)) n W u ([0, T]; (tf 1 ^))') D L 2 (0, T; H^Sl)), 
ceL°°(Q), 

given by 

fo(f,?)=fc°a,?) + eZ> 1 (M), 

(88) 

c(^)=c°(*.?) + ec 1 (M). 
In particular, the restrictions of the solution to each layer have the properties 
k G W u ([0, J]; (ff 1 ^,-))') n L°°(0, T;// 1 ^;)) n C([0, r];L 2 (£2,-)) n L 2 (0, T\H 2 (Qi)), 

ct e w 1 - 00 ^,-) n w^ao.rhff 1 ^,-)) n c^to.rjjL 2 ^,)) n w u ([0,r];ff 2 (£2,-)). 



4. The control problem 

Using the notation (60) we can rewrite the control problem (P) as 

Min(^ f ((c 0 (t,^ + sc l (t^)-c at ) + ) 2 d^dt+^ [ ((b°(t^)-b c ) + ) 2 d^dt) (P) 

b 0 eU\2J Q 2 Jq ) 

subject to Equations (53)-(59) (equivalently Equations (71) and (72)) and (84), where 

t/={/7oeL°°(^); b m <b 0 <b M }. 

This means that b° is the solution to the ^-approximation, Equations (53)-(59), with the initial 
datum bo and c is given by Equation (84). Here, b m — (&,,m)!=i,...,n and bu = (fei,M)(=i,...,«, and 
we note that U — Y\j=i Ui- 

We recall that c°(t, §) = cq and is constant in all layers. Because in (P) only the e°-order 
approximation for b and only the e 1 -order approximation for c are involved, for the writing 
simplicity we shall indicate them by b and c, without superscripts. 



Journal of Biological Dynamics 



101 



4.1. Existence of the optimal control 

Theorem 4. 1 Problem (P) has at least a solution. 

Proof Since J (bo) > 0, its infimum exists and let us denote it by 

d= inf ( l - f (( Co + ec ( f ,§)- Ccrt )+) 2 d§df+^ f {{b{t^)-b c )+) 2 d^ dt) >0. 

Let us take a minimizing sequence (bo)/t>i> b\ G U, such that 

d < Q jT (( Co + ec *( f , £) - Ccrt )+) 2 d| dr + | ((**('. §) " ^) + ) 2 d| df) < d + p (89) 

where b k = (b k ),-=i,...,„ is the solution to Equations (53)-(59) with initial datum b% and c* is given 
by Equation (84) corresponding to b k . We stress that the confusion between the sequence (b k ) k >i 
and the e A -order approximation should be avoided. 

Since b\ G U we can select a subsequence (denoted still by k) such that 

b k 0 -> fog weak-star in L°°(£2), as £ oo. 

Then, the Cauchy problem (71) and (72) 

dM t 

— (0 + Afe* (t) = 0 a.e. f G (0, T) , 
dt 

b k (0) = b k 

has, by Theorem 3.1, a unique solution b k , satisfying Equations (80)-(82) and (85). Therefore, 
selecting a subsequence we deduce that 

b k -> b* weakly in L 2 (0,T;V), as £ -> oo, 

dfo* d£* , , 

> weakly in L 2 (0,T;L 2 (Q)), as k ->• oo, 

df df 

which implies by Lions' theorem (since V is compact in L 2 (£2)) that 

b k ->• strongly in L 2 (0, 7;L 2 (^)), as £ ->• oo. 

Then, M — > b* a.e. on <2 and since ^ is continuous it follows that (p(b k (r , §), c o) - ► ^(&* > I). Co) 
a.e. on g, implying that (c 0 + sc k — c crt ) + — > (c 0 + sc* — c crt ) + a.e. in Q, as k — > oo. Moreover, 
^ is bounded by (is ), so by the Lebesgue dominated convergence theorem, we get by Equations (84) 
that 

(c 0 + ec* - c crt ) + -> (c 0 + ec* - c crt ) + strongly in L 1 (Q), as £ -> oo. 

But, on the other hand, by Equation (89), the sequence (co + ec k — c crt )^ is bounded in L 2 (Q), 
so that it converges weakly on a subsequence to a limit in L 2 (Q), and in conclusion, by the limit 
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uniqueness 

(cq + sc k — c crt ) + — > (co + ec* — c crt ) + weakly in L 2 (Q), as k — > oo. 

Since (fM(b k ))k lies in a bounded subset of L 2 (Q) and /x(fe*) — > n(b*) a.e. in g, it follows that 

ti(b k ) -> /z(Z>*) weakly inL 2 (g), as £ — >■ oo. 

We also deduce by the Ascoli-Arzela theorem that 

b k {t) -* b*(t) strongly in L 2 (Q), uniformly in t e [0, T], as k -> oo, 

which implies that M(0) = &q -> fog as k -> oo. Moreover, is the solution to Equations (71) 
and (72) because, on the basis of the above convergencies, we can pass to the limit as k — > oo in 
the weak form 



IQ 

getting 

'db 



db k \ 

— ir+ D(x)Vb k ■ Vir + ti{b k )f d§ df = 0 Vf e L 2 (0, T; V) 
at I 



LC 



ir + D(x)Vb* • Vf + ix(b*)f d£ df = 0 <= Z/(0, T; V). 
at ) 



Finally, using the weakly lower semicontinuity property of the functions in Equation (89), we 
can pass to the limit and deduce 



d< 1 - I ((c 0 + sc*(t,^-c cn ) + ) 2 d^dt+^ f ((b*(t^)-b c ) + ) 2 d^dt 
1 JQ 1 JQ 

< limmf Q j ((c 0 + sc k (t, £) - c crt )+) 2 d| dt + °- J ((b k (t, §) - b c )+) 2 d| dr) 



'G z -/G 

which shows that J{b* Q ) — d, hence (P) has at least a solution. ■ 
4.2. The system in variations and the dual system 

Let (b*Q , b* , c* ) be an optimal pair in (P) . The function />* is the solution to Equations (71) and (72) 
with the initial datum b* Q and c* is given by Equation (84) corresponding to b* , 



c*(t,i;)= [ <p(b*(s,^,c 0 )ds. 
Jo 

Let X > 0 and introduce b$ — b\* + Xw, where 



w — v — b* Q , v G U. 

The corresponding solution to Equations (71) and (72) with the initial datum b^ is denoted by b x , 
and we define 

bHt,H)-b*{t,H) 

Y(t,£) — hm . 



Journal of Biological Dynamics 



103 



By a few computation, it can be shown that Y solves the system in variations (written for F,- in 
each layer) 



dt 



DDjAYi + a*(t,t;)Yi = 0, (*,?)€&, i = l,n, 
y ( (0,f) = w h % e Q h i=hn, 
= 0, (f.fOeE,-, 



- 97,- - 3y i+ A 

9x 3x / 



-DD 



Yi{t, Xi ,H') = Yi+i(t,Xi,?), (MO e S,-, i = l,n - 1, 
9Fi 



dx 

djn 

dx 



0, e S 0 . 



x=x 0 



= 0, (a')eE», 



(90) 
(91) 
(92) 
(93) 
(94) 

(95) 
(96) 



VY t ■ v = 0, on E, lat , i = 1, 
where vv ; = v,- - b* 0 , with v,- G [/,■ and a* : (0, T) x £2, — > M, 

flT(M) = o*.-)* l (6r(M),c W ). 

We see that the system in variations can be treated in the same manner as the system for the e 1 -order 
approximation (being much simpler than it), and so it follows that it has a unique solution 

Y t e W 1 ' 2 ([0,r];(// 1 (fi / ))')nL 2 (0,r;// 1 (^«))nC([0,r];L 2 (^)). 
We introduce now the system for the dual variable p, written for each layer, 
dpi 



dt 



+ DDiApi-a*(t^)Pi 

= s^(b*(t,^,c 0 ) J (co + ec *( ff ,§) - Ccrt )+ ds + o{b*- b s y 



9ri 

(t^)eQi, i=l,n, 
Pi(r,f) = 0, |e£2i, i = T^, 



9x 9x 



= 0, (a'JeS, 



Pi(t,Xi,?) = p i+ i(t,Xi,$'), (MO G i=\,n-\, 
dpi 



-DDi 
DD n 



dx 

dPn 



= 0, (MOeSo, 



9x 



= 0, (a')eS„ 

! 

V/7; • v = 0, on S,! at , i = 1,.. .,«. 
It is obvious that it has a unique solution 



(97) 
(98) 

(99) 
(100) 

(101) 
(102) 

(103) 
(104) 



Pi g w u ([o,ri;(ff 1 (n i ))0nL 2 (o,r;ff 1 (£2 i ))nc([o,ri;z, 2 (n I -)). 
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4.3. Optimality conditions 

Proposition 4.2 Let (b^, b*) be optimal in (P). Then, the optimality conditions read 
b* 0 = b itm , on{£e£2 i ;/> j (0,t)>0}, 

bt 0 e (fo,-, m , b iM ) on G £2,-; Pi (0, £) = 0}, (105) 
b* 0 = b iM , on{f e£2,-;p i (0,?)<0}. 

Proo/ We write that (fog, &*) is optimal, that is, 

^ / ( 6 ') for an y 0 e t 7 - 

We take Q — b x and have 

^ f ((c 0 + £c*(a)-c crt )+) 2 d£df+^ f ((fo*(f,?)-fo c )+) 2 d^dr 
I Jo 2- Jo 



< I /" ((co+«r*(a)-c crt ) + ) 2 dfdf+^ /" §)-fc c )+) 2 d$df; 
2 Jo 2 Jo 



whence, passing to the right-hand side, performing some computations, dividing by X and letting 
X to go to 0, we get 

[ (c 0 + sc*(t,i;)-c cn ) + e [ <p ri (b*(a,i;),co)Y(s,$)ds + a [ (b*(t,$) - b s )+Y d£ dt > 0, 
Jq Jo Jq 

where <p n generically indicates d<Pi/dr\ in each layer. This can still be written, after changing the 
order of integrations, as 

f e<p n (b*(t,l;),co)Y(t,%) f (c 0 +ec*(s,i;)-c cn ) + ds + cT f (b*(t,$)-b s ) + Ydi;dt>0. 
Jq Jt Jq 

(106) 

Then, we test Equation (90) by p t and integrate over (0, T), using the initial data and boundary 
conditions. We obtain writing globally 



J Je^ ri (b*(t, £), c 0 ) jf (c 0 + sc*(s, £) - c crt ) + ds + a(b*(t, §) - b s )+ J 7 d| df 
- f *>(<M 

./ft 



| (/;"(/. £). c 0 ) / ^ 

,tM0,f)df=0, (107) 



where we recall that w(f) = v(|) — fog(t)- Combining Equations (106) and (107) we finally 
deduce 



[ p(0,$)(vG)-*5(£))d!>o, 
Jn 



which can be still written 

-P(P,& e dlu&t) = N a (PS). (108) 

In this relation, dluib^) is the subdifferential at b^ of the indicator set of U and NuQ}^) is the 
normal cone to t/ at b$. In conclusion, we obtain relations (105) as claimed. 
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We recall that b m and b M represent a sequence of minimum and maximum values assigned 
for each layer i, Equivalently, the previous relation — p(Q, f) e Nt/^o) means that — pt(0, §) e 
^V(7 f (^*o) an< i so tne relations given by relations (105) follow. ■ 

We observe that the minimization problem admits at least a solution which can be realized not 
necessarily at bo — b m , but in a way expressed by relations (105). We conclude that the controller 
/?* (= b* Q ) may take either the value £>, m or £>, m on the subsets of f2, , where pi (0, £ ) has positive or 
negative values. In some layers, b* 0 may reach also the value 0 (if e.g./?,(0, f ) > 0 in ft;), because 
£>,; m was allowed to be nonnegative. This means that in some situations it might be possible not 
to place bacteria in some layers at the initial time. Finally, we assert that it is difficult to say if the 
system (98)-(104) can have a null solution p,(0, £) on a subset of £2,-, but this analysis is beyond 
the scope of this paper. Therefore, we cannot exclude that in some circumstances it may happen 
that the minimum be reached for whatever bjo e (£> !;II1 , £>,,m)- 



5. Numerical results 

Some graphics revealing the feature of the process with respect to the change of bo are presented. 
The simulations are made with Comsol Multiphysics (FLN License 1025226) for the system 
(6)— (17) in a ID domain with the same values D, 8, /, cq in all layers (in order to compare how 
results change for various bo) with the following dimensionless data: 

(xo, Xn ) = (0, 10), T = 5, D=l, /=1, 3 = <p = e = 0.01, K = x = 5, 

3 = 1, c 0 = 10, 



D = 1, /(&, c) - -6, c) = x 



(1 +c) 2 ' 



c) = —be. 



The figures represent the values of b and c along time, drawn at x fixed (0, 2, 4, 6, 8, 10) 
corresponding to two values of the initial datum bo- In Figure 1, these graphics are plotted for bo, 
the step function 

20, xe[0,3), 
bo{x) =10, x e [3,6), (109) 
5, jc e [6, 10], 



Solution i at x fixed, for chi=5, bo step function 



Solution c at x fixed, for chi=5, f>o step function 




Figure 1. Solutions b(t, §) (left) andc(?,§) (right) for 6 0 W given by (109). 
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while in Figure 2 the graphics are displayed for 



A)inv 



(X) 



0, xe[0,3), 
20, .re [3,6), 
5, * e [6, 10] 



(HO) 



In Figure 1 (right), we see that the values of c increase from x — 0 to x — 10. The medium 
cleaning is more efficient beginning with the first layers. In Figure 2 (right), the medium cleaning 
happens in a different way, being less efficient in the first layers (x e [0, 2]) where the values of 
c remain high, but more efficient in the medium layers (x e [4, 6]). 



6. Conclusion 

We address a control problem related to a chemotactic motion of a bacteria towards a pollutant 
chemoattractant in a stratified medium. The control problem is focused on the reduction of the 
chemoattractant concentration, by acting upon the initial distribution of the bacteria population, 
bo. We prove that the control problem has at least a solution and provide the structure of the 
controller bo in each layer. 
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